IEEE TRANSACTIONS ON SIGNAL PROCESSING (TO APPEAR) 



1 



Distributed Optimal Beamformers for Cognitive 
Radios Robust to Channel Uncertainties 

Yu Zhang, Student Member, IEEE, Emiliano Dall'Anese, Member, IEEE, 
and Georgios B. Giannakis, Fellow, IEEE 



o- 

CO. 



> 

o 

oo 



OV 

o: 



X 



Abstract — Through spatial multiplexing and diversity, multi- 
input multi-output (MIMO) cognitive radio (CR) networks can 
markedly increase transmission rates and reliability, while con- 
trolling the interference inflicted to peer nodes and primary users 
(PUs) via beamforming. The present paper optimizes the design of 
transmit- and receive-beamformers for ad hoc CR networks when 
CR-to-CR channels are known, but CR-to-PU channels cannot be 
estimated accurately. Capitalizing on a norm-bounded channel 
uncertainty model, the optimal beamforming design is formulated 
to minimize the overall mean-square error (MSE) from all data 
streams, while enforcing protection of the PU system when the 
CR-to-PU channels are uncertain. Even though the resultant 
optimization problem is non-convex, algorithms with provable 
convergence to stationary points are developed by resorting to 
block coordinate ascent iterations, along with suitable convex 
approximation techniques. Enticingly, the novel schemes also lend 
themselves naturally to distributed implementations. Numerical 
tests are reported to corroborate the analytical findings. 

Index Terms — MIMO wireless networks, cognitive radios, 
beamforming, channel uncertainty, robust optimization, dis- 
tributed algorithms. 



I. Introduction 

Cognitive radio (CR) is recognized as a disruptive technol- 
ogy with great potential to enhance spectrum efficiency. From 
the envisioned CR-driven applications, particularly promising 
is the hierarchical spectrum sharing [1|, where CRs oppor- 
tunistically re-use frequency bands licensed to primary users 
(PUs) whenever spectrum vacancies are detected in the time 
and space dimensions. Key enablers of a seamless coexistence 
of CR with PU systems are reliable sensing of the licensed 
spectrum (2), [3), and judicious control of the interference 
that CRs inflict to PUs [1 1. In this paper, attention is focused 
on the latter aspect. 

Recently, underlay multi-input multi-output (MIMO) CR 
networks have attracted considerable attention thanks to their 
ability to mitigate both self- and PU-inflicted interference via 
beamforming, while leveraging spatial multiplexing and diver- 
sity to considerably increase transmission rates and reliability. 
On the other hand, wireless transceiver optimization has been 
extensively studied in the non-CR setup under different design 
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criteria (4), (5], and when either perfect or imperfect channel 
knowledge is available; see e.g., J6), (7), and references 
therein. In general, when network-wide performance criteria 
such as weighted sum-rate and sum mean-square error (MSE) 
are utilized, optimal beamforming is deemed challenging 
because the resultant optimization problems are typically non- 
convex. Thus, solvers assuring even first-order Karush-Kuhn- 
Tucker (KKT) optimality are appreciated in this context |@)- 

In the CR setup, the beamforming design problem is ex- 
acerbated by the presence of interference constraints (TJ. In 
fact, while initial efforts in designing beamformers under 
PU interference constraints were made under the premise 
of perfect knowledge of the cognitive-to-primary propagation 
channels [8|-[11|, it has been recognized that obtaining ac- 
curate estimates of the CR-to-PU channels is challenging or 
even impossible. This is primarily due to the lack of full 
CR-PU cooperation (T), but also to estimation errors and 
frequency offsets between reciprocal channels when CR-to-PU 
channel estimation is attempted. It is therefore of paramount 
importance to take the underlying channel uncertainties into 
account, and develop prudent beamforming schemes that en- 
sure protection of the licensed users. 

Based on CR-to-PU channel statistics, probabilistic inter- 
ference constraints were employed in lfl2l for single-antenna 
CR links. Assuming imperfect knowledge of the CR-to- 
PU channel, the beamforming design in a multiuser multi- 
input single-output (MISO) CR system sharing resources with 
single-antenna PUs was considered in [13|; see also [ 14 1 for a 
downlink setup, where both CR and PU nodes have multiple 
antennas. The minimum CR signal-to-interference-plus-noise 
ratio (SINR) was maximized under a bounded norm constraint 
capturing uncertainty in the CR-to-PU links. Using the same 
uncertainty model, minimization of the overall MSE from 
all data streams in MIMO ad hoc CR networks was con- 
sidered in [ 1 5 1 . However, identical channel estimation errors 
for different CR-to-PU links were assumed. This assumption 
was bypassed in [16], where the mutual information was 
maximized instead. Finally, a distributed algorithm based on 
a game-theoretic approach was developed in ifTTl . 

The present paper considers an underlay MIMO ad hoc 
CR network sharing spectrum bands licensed to PUs, which 
are possibly equipped with multiple antennas as well. CR- 
to-CR channels are assumed known perfectly, but this is not 
the case for CR-to-PU channels. Capitalizing on a norm- 
bounded uncertainty model to capture inaccuracies of the CR- 
to-PU channel estimates, a beamforming problem is formu- 
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lated whereby CRs minimize the overall MSE, while limiting 
the interference inflicted to the PUs robustly. The resultant 
robust beamforming design confronts two major challenges: 
a) non-convexity of the total MSE cost function; and, b) the 
semi-infinite attribute of the robust interference constraint, 
which makes the optimization problem arduous to manage. 
To overcome the second hurdle, an equivalent re-formulation 
of the interference constraint as a linear matrix inequality 
(LMI) is derived by exploiting the S-Procedure fl8l . On the 
other hand, to cope with the inherent non-convexity, a cyclic 
block coordinate ascent approach |fl9l is adopted along with 
local convex approximation techniques. This yields an iterative 
solution of the semi-definite programs (SDPs) involved, and 
generates a convergent sequence of objective function values. 
Moreover, when the CR-to-CR channel matrices have full 
column rank, every limit point generated by the proposed 
method is guaranteed to be a stationary point of the original 
non-convex problem. However, CR links where the transmit- 
ter is equipped with a larger number of antennas than the 
receiver, or spatially correlated MIMO channels [20|, can 
lead to beamformers that are not necessarily optimal. For this 
reason, a proximal point-based regularization technique ED is 
also employed to guarantee convergence to optimal operating 
points, regardless of the channel rank and antenna configu- 
ration. Similar to H-fli], ED-El, E2-C1, perfect time 
synchronization is assumed at the symbol level. 

Interestingly, the schemes developed are suitable for dis- 
tributed operation, provided that relevant parameters are ex- 
changed among neighboring CRs. The algorithms can also be 
implemented in an on-line fashion which allows adaptation to 
(slow) time-varying propagation channels. In this case, CRs 
do not necessarily wait for the iterations to converge, but 
rather use the beamformer weights as and when they become 
available. This is in contrast to, e.g., flU, ||6) and 1 14l- lfl6l 
in the non-CR and CR cases, respectively, where the relevant 
problems are solved centrally and in a batch form. 

In the robust beamforming design, the interference power 
that can be tolerated by the PUs is initially assumed to be 
pre-partitioned in per-CR link portions, possibly according 
to quality-of-service (QoS) guidelines ifTTI . lUTl . However, 
extensions of the beamforming design are also provided when 
the PU interference limit is not divided a priori among CR 
links. In this case, primal decomposition techniques [19| are 
invoked to dynamically allocate the total interference among 
CRs. Compared to [15|, the proposed scheme accounts for 
different estimation inaccuracies in the CR-to-PU links. 

The remainder of the paper is organized as follows. Sec- 
tion [n] introduces the system model and outlines the proposed 
robust beamforming design problem. In Section [HI] the block 
coordinate ascent solver is developed, along with its proximal 
point-based alternative. Aggregate interference constraints are 
dealt with in Section [IV] and numerical results are reported 
in Section [V] Finally, concluding remarks are given in Sec- 
tion [VT] while proofs are deferred to the Appendix. 

Notation. Boldface lower (upper) case letters represent 
vectors (matrices); H nxn , H" xn , C nxn and K stand for spaces 
ofnxn Hermitian, n x n Hermitian positive semidefinite, 
n x n complex matrices, and real numbers, respectively; (-) T , 




Fig. 1. The system model for MIMO ad hoc CR networks. 

(•)*, and {-) u indicate transpose, complex conjugate, and 
conjugate transpose operations, respectively; Tr{-} denotes the 
trace operator, and vec(A) the vector formed by stacking the 
columns of A; ||a||2 and |A||_f represent the Euclidean norm 
of a and the Frobenius norm of A, respectively; A <g) B is 
the Kronecker product of A and B; Ijv is the NxN identity 
matrix; Finally, E{-} denotes the expectation operator, and 
K(-) stands for the real part of a complex number. 

II. System Model and Problem Formulation 

Consider a wireless MIMO CR network comprising K 
transmitter-receiver pairs {Uf,, Ul}. Let and Nf., k G 
JC := {1,2,..., K}, denote the number of antennas of the fc-th 
transmitter-receiver pair, as shown in Fig.Q] Further, let Sfc de- 
note the Mfc x 1 information symbol vector transmitted by L7| 
per time slot, with covariance matrix Ejsfc.s^} = In or- 
der to mitigate self-interference, transmitter Ul pre-multiplies 
Sk by a transmit-beamforming matrix Ffc £ C MfeXMfe ; that is, 
Ul actually transmits the x 1 symbol vector Xfc := FfcSfc. 
With Hfcj G C NkXM i denoting the Uj to U 7 k channel matrix, 
the Nk x 1 symbol received at UJ1 can be written as 

y k = H fc!fc x fe + ^2 H feJ x j + n fe (!) 

where G C Nk is the zero-mean complex Gaussian dis- 
tributed receiver noise, which is assumed independent of Sfc 
and {Hfc j}, with covariance matrix E{nfcn^} = <r^Ijv fe . 

Low-complexity receiver processing motivates the use of a 
linear filter matrix Wfc G C MkXNk at UL to recover S& as 

s fe := W fc y fe , k G K. (2) 

Using W fc at UL the MSE matrix E fc := 
E{(§fc — Sfc) (sk — Sk) }, which quantifies the reconstruction 
error, is given by [cf. ([U] 

E fc = W,A fc W« - WfeH^Ffe - F«H« fe W« + I Mk 

(3) 

where A k := J^f =1 TT^F.F^n^ + ^I^. Entry of 
Efc represents the MSE of the i-th data stream (i-th entry of 
Sfc) from Ul to UL and Tr{Efc} corresponds to the MSE of 

Sfe. 
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To complete the formulation, let Gfc £ i£ LxM k denote 
the channel between CR Ut and a PU receiver, possibly 
equipped with multiple (L) antennas^, and i max the maximum 
instantaneous interference that the PU can tolerate. As in 
e-g-, irm . 02], suppose that t max is pre-partitioned in per- 
CR transmitter portions {i™ ax }, possibly depending on QoS 
requirements of individual CR pairs. Then, the transmit- and 
receive-beamforming matrices minimizing the overall MSE 
can be obtained as (see also 0) 

K 

(PI) min VTr{E fc } (4a) 

s. t. Tr{F fc F«} < P r x , ke)C (4b) 
Tr{G fc F fc F«G«} < k £ K (4c) 

where p™ ax is the maximum transmit-power of Ufa. 
Remark 1 (Adopted performance metric). Among candidate 
performance metrics, the sum of MSEs from different data 
streams is adopted here, which has been widely employed in 
the beamforming literature; see e.g., (6), Q and references 
therein. The relationships between MSE, bit error rate (BER) 
and SINR have been thoroughly considered in 1221 . and further 
investigated in [6|. Specifically, it has been shown that an 
improvement in the total MSE naturally translates in a lower 
BER. Furthermore, the sum of MSEs facilitates derivation of 
optimal filters, and the equivalence between minimizing the 
weighted sum of MSEs and maximizing the weighted sum 
rate has been established in [4], |23|. 

Unfortunately, due to lack of explicit cooperation between 
PU and CR nodes, CR-to-PU channels {Gfc} are in general 
difficult to estimate accurately. As PU protection must be 
enforced though, it is important to take into account the CR- 
to-PU channel uncertainty, and guarantee that the interference 
power experienced by the PU receiver stays below a prescribed 
level for any possible (random) channel realization [12], |14|. 
Before developing a beamforming approach robust to inac- 
curacies associated with channel estimation, problem (PI) is 
conveniently re-formulated first in order to reduce the number 
of variables involved. 

A. Equivalent Optimization Problem 

For the sum-MSE cost in (f4ah . the optimum {W^ pt } will 
turn out to be expressible in closed form. To show this, note 
first that for fixed {F fc }, (PI) is convex in W fc , and {W° pt } 
can be obtained from the first-order optimality conditions. 
Express the Lagrangian function associated with (PI) as 

K K 

C (V, V)=J2 Tr{E fc } + E ^ (Tr{F fc F«} - p^) 

k=l k=l 
K 

+ 5> (Tr{G fe F fe F*G*} ~ %*) (5) 

fc=i 

where V := {F fc ,W fc }f =1 and V := {A fe ,i/ fe }f =1 collects 
the primal and dual variables, respectively. Then, by equating 

'A single PU receiver is considered throughout the paper. However, 
extension to multiple receiving PUs is straightforward; see also Remark 5. 



the complex gradient DC ("P,2?) /<9W£ to zero, matrix W^ pt 
is expressed as 

Wr = F«H« fe A-\ fce/C. (6) 

Clearly, the optimal set {W^ pt } does not depend on channels 
{Gfc}, but only on {H fcj }. 

Substituting {W^ pt } into (Bal l, and using the covariance 
Qfc := E{xfcX^} = FfcF^ as a matrix optimization variable, 
it follows that (PI) can be equivalently re- written as 

K 

(P2) max Vufc({Qfc}) (7a) 

s. t. Tr{Q fc } < p£», k e K (7b) 
Tr{GfcQ fc G«} < lT\ kelC (7c) 

where the per-CR link utility ({Qfc}) is given by 

u k ({Qfc}) := Tr{Hfc,fcQfcH« fc ( I I^-Qa 1 1 + Rfc^fc)" 1 } 

(8) 

with R fej fc := Yi&k H mQ* H m + ^Un,,- One remark is in 
order regarding (P2). 

Remark 2 (Conventional MIMO networks). Upon discarding 
the interference constraints d7cl >. the beamforming problems 
formulated in this paper along with their centralized and 
distributed solvers can be considered also for non-CR MIMO 
ad-hoc and cellular networks in downlink or uplink operation. 

Channels {Gfc} must be perfectly known in order to solve 
(P2). A robust version of (P2), which accounts for imperfect 
channel knowledge, is dealt with in the next section. 

B. Robust Interference Constraint 

In typical CR scenarios, CR-to-PU channels are challenging 
to estimate accurately. In fact, CR and PU nodes do not 
generally cooperate [1], thus rendering channel estimation 
challenging. To model estimation inaccuracies, consider ex- 
pressing the CR-to-PU channel matrix Gfc as 

Gfc = Gfc + AGfc , fce/C (9) 

where Gfc is the estimated channel, which is known at CR 
transmitter U%, and {AGfc} captures the underlying channel 
uncertainty fl4l . fl6l . Specifically, the error matrix AGfc is 
assumed to take values from the bounded set 

C7fc := {AGfc |Tr{ AGfc AG^} < e 2 k } , k e IC (10) 

where efc > specifies the radius of Qk, and thus reflects 
the degree of uncertainty associated with Gfc. The set in (fTOl i 
can be readily extended to the general ellipsoidal uncertainty 
model |[T8l Ch. 4]. Such an uncertainty model properly 
resembles the case where a time division duplex (TDD) 
strategy is adopted by the PU system, and CRs have prior 
knowledge of the PUs' pilot sequence(s). But even without 
training symbols, CR-to-PU channel estimates can be formed 
using the deterministic path loss coefficients, and the size of 
the uncertainty region can be deduced from fading channel 
statistics. Compared to [12], the norm-bounded uncertainty 
model leads to worst-case interference constraints that ensure 
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PU protection for any realization of the uncertain portion of 
the propagation channels. 

Based on (TlOb . a robust interference constraint can be 
written as 

Tr{(G fc + AG fc )Q fe (G fe + AG/} < if, 

VAG t e&, fce/C (11) 

and consequently, a robust counterpart of (P2) can be formu- 
lated as follows 



(P3) max > u k ({Qk}) 



K 



{Qfc>:0} 



fe=l 



(12a) 



s. t. Tr{Q fe } < P f, fce/C (12b) 
Tr{G fc Q fc G«} < if, V AG fc g Q k , fce/C. (12c) 

Clearly, once {Q^ pt } are found by solving (P3), the wanted 
{W^ pt } can be readily obtained via (|6), since Afc := 

12j=i H fc,iQj' H ^j -+°l l N k - However, Y. k u k ({Qk}) is non- 
convex in {Qfc}, and hence (P3) is hard to solve in general. 
Additionally, constraints ( fTTT i are not in a tractable form, which 
motivates their transformation. These issues are addressed in 
the next section. But first, two remarks are in order. 
Remark 3 (Uncertain MIMO channels). In an underlay hier- 
archical spectrum access setup, it is very challenging (if not 
impossible) for the CRs to obtain accurate estimates of the 
CR-to-PU channels. In fact, since the PUs hold the spectrum 
license, they have no incentive to feed back CR-to-PU channel 
estimates to the CR system [1). Hence, in lieu of explicit CR- 
PU cooperation, CRs have to resort to crude or blind estimates 
of their channels with PUs. On the other hand, sufficient time 
for training along with sophisticated estimation algorithms 
render the CR-to-CR channels easier to estimate. This explains 
why similar to relevant works IfTZl - lfTTll , CR-to-CR channels 
are assumed known, while CR-to-PU channels are taken as 
uncertain in CR-related optimization methods. Limited-rate 
channel state information that can become available e.g., with 
quantized CR-to-CR channels (6), 0, 1121 . can be considered 
in future research but goes beyond the scope of the present 
paper. 

Remark 4 (Radius of the uncertainty region). In practice, 
radius and shape of the uncertainty region have to be tailored 
to the specific channel estimation approach implemented at 
the CRs, and clearly depend on the second-order channel 
error statistics. For example, if AGt has zero mean and 
covariance matrix Sag,. = oil, where u\ depends on the 
receiver noise power, and the transmit-power of the PU (see, 
e.g. [25|), then the radius of the uncertainty region can be set 
to t\ = K^k&k' wnere £fc denotes the path loss coefficient, 
and k > a parameter that controls how strict the PU 
protection is. Alternatively, the model e\ — Kaf.\\Qk\\ 2 p can 
be utilized ifTTl . If Sag,. ^ oil, then the uncertainty region 



can be set to Q k := {AG fe |Tr{AG fc S^ fc AG^} < n) 
and the robust constraint (fTTT i can be modified accordingly. 
Similar models are also considered in [261. 



III. Distributed Robust CR Beamforming 

To cope with the non-convexity of the utility function dl2at . 
a block-coordinate ascent solver is developed in this sec- 



tion. Define first 

A(Qfc,Q-fc) := 



the sum of all but the 



with Q _ fc 



fc-th utility as 

{Qoli + *}■ 



Notice that Ufc(-) is concave and /&(•) is convex in Qfc; see 
Appendix [B] for a proof. Then, (P3) can be regarded as a 
difference of convex functions (d.c.) program, whenever only 
a single variable Qfc is optimized and Q_fc is kept fixed. 
This motivates the so-termed concave-convex procedure |27|, 
which belongs to the majorization-minimization class of al- 
gorithms [28 1, to solve problem (P3) through a sequence of 
convex problems, one per matrix variable Qfc. Specifically, the 
idea is to linearize the convex function /&(•) around a feasible 
point Qfc, and thus to (locally) approximate the objective ( 112at 
as (see also [5| and [10|) 



K 

E 

k=\ 



u k ({Q fe }) = ma.. ({Q fc }) + / fc (Qfc, Q- fe ) 
/ fe (Q fc ,Q-fc)+Tr{D«(Q A 



~ u k ({Qfc}) 
where 



Qfe)} (13) 



V Qfc / fc (Q fc ,Q_ fc ) 



df k 



9Qt 



(14) 



Q fc =Q fc 

Therefore, for fixed Q_fc, matrix Qfc can be obtained by 
solving the following sub-problem 



(P4) max U fc(Qfc,Q_fc)+Tr{D«Qfc} 

s. t. Qfc y o 

Tr{Qfc}< P r 
Tr{GfcQfcG«} < il 

where (see Appendix lAl 



V AGfc e Qk 



D * :=-E H ,W V i B 7 lH 



Qfc=Q fe 



K 



B i : E " ( 2 " •'< + vfixy, 

i=l 



(15a) 

(15b) 
(15c) 
(15d) 



(16) 



(17) 



(18) 



At each iteration n = 1,2,..., the block coordinate ascent 
solver amounts to updating the covariance matrices {Qfc} in 
a round robin fashion via (P4), where the solution obtained at 
the (n — l)-st iteration are exploited to compute the complex 
gradient ( TTol ). The term TrjD^Qfc} discourages a "selfish" 
behavior of the fc-th CR-to-CR link, which would otherwise 
try to simply minimize its own MSE, as in the game-theoretic 
formulations of IfTTl and ifTTll . In the next subsection, the 
robust interference constraint will be translated to a tractable 
form, and (P4) will be re-stated accordingly. 

A. Equivalent robust interference constraint 

Constraint (1 1 5db renders (P4) a semi-infinite program 
(cf. ||29l Ch. 3]). An equivalent constraint in linear matrix 
inequality (LMI) form will be derived next, thus turning (P4) 
into an equivalent semi-definite program (SDP), which can be 
efficiently solved in polynomial time by standard interior point 
methods. To this end, the following lemma is needed. 
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Lemma 1. (S-Procedure H8\ p. 655]) Consider A,D £ 
H™ xn ,b € C,c.e € K, fl«rf assume the interior condition 
holds, i.e., there exists an x satisfying x^Dx < e. 77zen, ?/ze 
inequality 



x H Ax + 25R(b^x) + c > 0, V x H Dx < e 
/zo/cfc if and only if there exists 6 > 5mc/i f/iaf 

b 

c — e£ 



(9D + A 
b* 



y o. 



(19) 



(20) 



Using Lemma [T] the robust constraint d!5db can be equiv- 
alently reformulated as follows. 

Proposition 1. There exists 9 k > 0, so that the robust 
interference constraint \15d\ is equivalent to the following LMI 



~>k i -LxM h 



{II ® Qfc) 



-vec(Q«G«)« 



-vec(Q«G«) 



, max 



-Tr{G fe Q fc G«} 



y o. 



(21) 



Proof: Using the properties of the trace operator 
Tr(Z w AZ) = vec(Z) w (I ® A)vec(Z) and Tr(B^Z) = 
vec(B) w vec(Z), constraint jl5dj can be re-written as 

- g« (I L ® Q fc ) gfc - 25R (vec(Q«G«) w gfe 



Tr{G fe QfcG^} > 0, V ||gfc|| 2 < efc (22) 

where gfc := vec(AG^). Then, applying Lemma Q] to ( l22l 
yields readily (f2Tb . ■ 

Proposition 2. Problem (P4) can be equivalently re-written 
as the following SDP form: 



(P5) min Tr{T}-Tr{D*Q fc } 

T,e fc >o 



s.t. Tr{Qfc}< P r 
Hfc,fcQfcH?^ fc + Rfc_fc R 



-vec(Q«Gj 



l fc,fc 

K fc,fc 

(II ® Qfc) 

Hr>HYH 



1/2 

T 



!- 



-vec(Q«G«) 



e|&fc 



-Tr{G fc Qfc.G«} 



(23a) 

(23b) 
(23c) 

y o. 



R 



fr.fr 



fW/: First, note that [cf. ©] 

Uk (Qfc, Q-fc) = Tr |ljv fc — Rfc.fc (Hfc ; fcQfcHfc^ fe 
Thus, (P4) is equivalent to 

™ n Tr { R ^fc (Hfc.fcQfcH^ + Rfc.fc)- 1 R^ 2 } 

-TrjD^Qfc} 

s. t. dBbJi - <fl5dl 



Then, an auxiliary matrix variable Y is introduced such 



that Y y R 



1/2 / 

k,k y 



H 



fc.fc 



QfeH^fc 



R 



fc.fc 



R 



1/2 
fe,*' 



which can 



be equivalently recast as (|23c| > by using the Schur comple- 
ment 1 1 8 1 . Combining the LMI form of the robust interference 
constraint (fJTJ, the formulation of (P5) follows immediately. 



Algorithm 1 Centralized robust sum-MSE minimization 
l: Collect all channel matrices {Hj fc}, and noise powers 

2: Collect all CR-to-PU channel matrices {Gfc}, and confi- 
dence intervals {efc} 
3: Initialize Q fc 0) = 0,V k G JC 
4: repeat (n = 1, 2, . . .) 
5: for fe = 1,2,..., K do 
6: Compute D fc n) via (O 

7: Update Qfc by solving (P5) [(P6) for the proximal 

point-based method] 
8: end for 

9: until U (Q (,l) ) - U (Q ( "- 1} ) < v 
10: Calculate {W° pt } via © 

11: Broadcast optimal transmit- and receive-beamformers 



Problem (P3) can be solved in a centralized fashion upon 
collecting CR-to-CR channels {H,-^}, CR-to-PU estimated 
channels {Gfc}, and confidence intervals {efc} at a CR fusion 
center. The optimal transmit-covariance matrices can be found 
at the fusion center by solving (P5), and sent back to all 
CRs. This centralized scheme is tabulated as Algorithm Q] 
where Qfc denotes the transmit-covariance matrix of CR 
U~l at iteration n of the block coordinate ascent algorithm; 
Q( n ) := ^Qi"^, . . . , Q^J represents the set of transmit- 
covariance matrices at iteration n; U(-) is the objective func- 
tion ( 1 1 2ab . A simple stopping criterion for terminating the 
iterations is U (Q^) —U (Q ( " -1 )) < v, where v > denotes 
a preselected threshold. 

To alleviate the high communication cost associated with 
the centralized setup, and ensure scalability with regards to 
network size and enhanced robustness to fusion center failure, 
a distributed optimization algorithm is generally desirable. 
It can be noticed that the proposed coordinate ascent ap- 
proach lends itself to a distributed optimization procedure 
that can be implemented in an on-line fashion. Specifically, 
each CR Uj. can update locally Qfc via (P5) based on a 
measurement of the interference R fc ' ^ flOl, and the following 
(23d) information necessary to compute the complex gradient ( fl6l i: 



i) its covariance matrix Q 
,(«) 



(n-1) 



obtained at the previous 

(n)l 



iteration; ii) matrices {B^ '} and {Vj ; } obtained from the 
neighboring CR links via local message passing. Furthermore, 
it is clear that the terms in (fTST l corresponding to CRs located 
far away from CR U\. are negligible due to the path loss 
effect in channel {H-j,k}jjtk\ hence, summation in ( fl6l l is 
only limited to the interfering CRs, and consequently, matrices 
{Bj } and {V^} need to be exchanged only locally. The 
overall distributed scheme is tabulated as Algorithm |2] The 
on-line implementation of the iterative optimization allows 
tracking of slow variations of the channel matrices; in this 
case, cross-channels {Hj.fc} in Algorithm [2] need to be re- 
acquired whenever a change is detected. Finally, notice that 
instead of updating the transmit-covariances in a Gauss-Seidel 
fashion, Jacobi iterations or asynchronous schemes |19| can 
be alternatively employed. 
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Algorithm 2 Distributed on-line robust sum-MSE minimiza- 
tion 



1: Initialize Q { k ' = 0,V k £ K, 
2: repeat (n = 1, 2, . . .) 
3: for k = 1,2,..., K do 

4: Uj. acquires Hj,fc from its neighboring UJ 

5: transmit {Bj, , } to neighboring nodes 

6: receive {B^™-* , }j^k from neighboring nodes 

7: Compute Dj^ via © 

8: Measure 

9: Update Qt by solving (P5) [(P6) for the proximal 

point-based method] 
10: Update W<f ' via © 

11: Transmit and receive signals using Q ( k n) and W ( k n) 

12: end for 

13: until u k (Q(«)) - u fc (Q^ 1 - 1 )) < u, V k e K 



Remark 5 {Multiple PU receivers). For ease of exposition, the 
formulated robust optimization problems consider a single PU 
receiver. Clearly, in case of TVpu > 1 receiving PU devices, or 
when a grid of /Vpu potential PU locations is obtained from the 
sensing phase [30], a robust interference constraint for each of 
the KNp V CR-to-PU links must be included in (P3). As for 
(P5), it is still an SDP, but with N PU LMI constraints (123dl >. 
and one additional optimization variable (Ok) per PU receiver. 
Remark 6 (Network synchronization). Similar to Bl- lfTTl . 
[13 1 — [ 17], I22I- I1241 . time synchronization is assumed to have 
been acquired. In practice, accurate time synchronization 
among the CR transmitters can be attained (and maintained 
during operation) using e.g., pairwise broadcast synchroniza- 
tion protocols H311 . consensus-based methods 11321 . or mutual 
network synchronization approaches [33|. To this end, CRs 
have to exchange synchronization beacons on a regular basis; 
clearly, the number of time slots occupied by the transmission 
of these beacons depends on the particular algorithm imple- 
mented, the CR network size, and the targeted synchronization 
accuracy. For example, the algorithm in ||3T1 entails two 
message exchanges per transmitter pairs, while the message- 
passing overhead of consensus-based methods generally de- 
pends on the wanted synchronization accuracy 11321 . Since 
the CR network operates in an underlay setup, this additional 
message passing can be performed over the primary chan- 
nels). Alternatively, a CR control channel can be employed 
to avoid possible synchronization errors due to the interference 
inflicted by the active PU transmitters. Analyzing the effect of 
mistiming constitutes an interesting research direction, but it 
goes beyond the scope and page limit of this paper. 

B. Convergence 

Since the original optimization problem (P3) is non-convex, 
convergence of the block coordinate ascent with local convex 
approximation has to be analytically established. To this end, 
recall that (P4) and (P5) are equivalent; thus, convergence 
can be asserted by supposing that (P4) is solved per Gauss- 
Seidel iteration instead of (P5). The following lemma (proved 



in Appendix IB1 is needed first. 

Lemma 2. For each k £ K,, the feasible set of problem (P4), 
namely Qk '■= {Qfe|Qfc £ (| 1 5fc[) — (|15<ip }, is convex. The 
real-valued function /fc(Qfc,Q— ft) is convex in over the 
feasible set Qk, when the set Q fc := {Qj;.? =/= k} is fixed. 

Based on Lemma [2] convergence of the block coordinate 
ascent algorithm is established next. 

Proposition 3. The sequence of objective function values d!2ab 
obtained by the coordinate ascent algorithm with concave- 
convex procedure converges. 

Proof: It suffices to show that the sequence of objective 
values d!2ab is monotonically non-decreasing. Since the ob- 
jective function value is bounded from above, the function 
value sequence must be convergent by invoking the mono- 
tone convergence theorem. Letting Uk{-) denote the objective 
function (I15ab , which is the concave surrogate of U(-) as the 
original objective ( 112al l. consider 

Qi n) := 

argmax U k ((fc; Q<"\ . . . , Q^, Q^, . . . , Q%~ l A 
Qk&Qk v 7 

(24) 

where n stands for the iteration index. Furthermore, define 

z („) ;= (Q (n+l) j ; Q (n+1) ) Q W 5 Q (n) }) (25) 

Q(n) ;= (Q („+l) j ; Q (n+ 1)j Q W 5 Q (n) } _ (2fi) 
Then, for all k £ K,, it holds that 

U (Z k n) ) = u k (Ct[ n+1 \ QW) + f h (Q k n+1 \ QW) (27a) 

>« fc (Q^Q5!)+A(Qi n) ,Qa) 

+ Tr{D«(Qi" +1) -Q( n) )} (27b) 

>u k (QW Q<S)+/*(Qir } ,QL B 2) 

+ Tr{Dj l (QW-QW)} (27c) 

= U (Z^) (27d) 

where (I27bb follows from the convexity of /fc(-) established in 
LemmalU (I27ct holds because Q[" +1 ' is the optimal solution 
of (P4) for fixed Q^. 

To complete the proof, it suffices to show that U (Q^™ +1 ^) 
is monotonically non-decreasing, namely that 

u (q ( " +1) ) > u (z£! x ) > . . . > u (zi n) ) > u (q (?1) ) 

(28) 
■ 

Interestingly, by inspecting the structure of {Hk,k, k £ /C}, 
it is also possible to show that every limit point generated by 
the coordinate ascent algorithm with local convex approxima- 
tion satisfies the first-order optimality conditions. Conditions 
on {Hfc fc,fc £ /C} that guarantee stationarity of the limit 
points are provided next. First, it is useful to establish strict 
concavity of the objective (|15at in the following lemma proved 
in Appendix ICl 



ZHANG et at: DISTRIBUTED OPTIMAL BEAMFORMERS FOR COGNITIVE RADIOS ROBUST TO CHANNEL UNCERTAINTIES 



7 



Lemma 3. If the channel matrices {H^^fe G /C} of the CR 
links {lit ~ > U[} have full column rank, then the objective 
function dl5at is strictly concave in 

We are now ready to establish stationarity of the limit points. 

Theorem 1. If matrices {H^, k 6 /C} have full column rank, 
then every limit point of the coordinate ascent algorithm with 
concave-convex procedure is a stationary point of (P3). 

Proof: The proof of Theorem Q] relies on the basic con- 
vergence claim of the block coordinate descent method in ||29l 
Ch. 2] and [5|. What must be shown is that every limit point 
of the algorithm satisfies the first-order optimality conditions 
over the Cartesian product of the closed convex sets. Let 
Q := (Qi, . . . , Qk) be a limit point of the sequence {Q^"-*}, 
and {Q^ nj ^\j = 1,2,...} a subsequence that converges to 
Q. First, we will show that lim Q^" J+1 - ) = Q\. Argue by 

contradiction, i.e., assume that {Q*" J+1) - Q^} does not 

converge to zero. Define 7^) := ||Q^" J+1) - QP } ||f- By 

possibly restricting to a subsequence of {nj}, it follows that 

there exists some 7 > such that 7 < 7"' for all j. 

Let S[ nj) := (Q^" J+1) - Q£ nj) )/7 (Ty) - Thus ' we have that 
Q K+i) = Q („,) +7 K) S K) and \\ S to)y = 1. Because 

(n ■ ) 

S\ belongs to a compact set, it can be assumed convergent 
to a limit point Si along with a subsequence of {nj}. 

Since it holds that < £7 < 7 (n ^ for all e £ [0, 1], the point 
Qi + tjSi lies on the segment connecting two feasible 
points Q^ j) and Q^ J+1) . Thus, Q^ j) + eyS^ is also 
feasible due to the convexity of Qi [cf. Lemma [2). Moreover, 
concavity of U\{-\ Q,_( ) implies that U\ is monotonically 

(n ■ ) 

non-decreasing in the interval connecting point Q\ to 
Q( n j+!) over tne set q x fj enC g ; \i readily follows that 

£i(Q^ +1) ; QL"/ ) > Ui(Q { i l3) + Q^) 

> WiCQ^jQ^). (29) 

Note that U\{-) is a tight lower bound of U(-) at each 
current feasible point. Also, from <|28), Ui(Q[ nj+1) ; Q^') 
is guaranteed to converge to Wi(Q) as j — > 00. Thus, upon 
taking the limit as j — > 00 in ( |29l , it follows that 

M 1 (Qi + e7Si;Q_ 1 )=W 1 (Q), V e £ [0, 1] . (30) 

However, since 7S1 7^ 0, (f30b contradicts the unique maxi- 
mum condition implied by the strict concavity of Ui(-\ •) in 
Qi [cf. Lemma [3). Therefore, Q^' lj+1 ^ converges to Qi as 
well. 

Consider now checking the optimality condition for Qi. 
Since Q^ n3+1 ' is the local (and also global) maximum of 
Wi(-;Q_i j) ), we have that 

SRjTrjv^! (q^' +1) ;Q^ 3 ) W (Qi - Q^ +1) ) }| < 0, 

VQiGQi (31) 

where ViUi(-) denotes the gradient of U%(-) with respect 
to Qi. Taking the limit as j 00, and using the fact that 



ViWi(Q) = ViWi(Q), it is easy to show that 

3?{Tr{ViW(Q) w (Qi-Qi)}} < 0, V Qi e Q X - (32) 
Using similar arguments, it holds that 

sR{Tr{ViW(Q) K (Qi-Qi)}} <0, 

VQ, £ Qi, i = l,2,...,X (33) 

which establishes the stationarity of Q and completes the 
proof. ■ 

C. Proximal point-based robust algorithm 

The full column rank requirement can be quite restrictive 
in practice; e.g., if M k > N k for at least one CR link, or 
in the presence of spatially correlated MIMO channels 1201 . 
Furthermore, computing the rank of channel matrices increases 
the computational burden to an extent that may not be afford- 
able by the CRs. In this section, an alternative approach based 
on proximal-point regularization II2TI is pursued to ensure 
convergence, without requiring restrictions on the antenna 
configuration and the channel rank. 

The idea consists in penalizing the objective of (P4) using 
a quadratic regularization term ^HQfc — Qjj™ \\%, with a 
given sequence of numbers 77 > 0. Then, (P5) is modified as 

(P6) min Tr{T} -Tr{D^Q fc } + — Tr{Y} (34a) 

T Y 

„ , *M k Qk - Qi" x) , n 

[Q fc -Qr i} y j-° 

(34b) 

d23bl. (I23cl (l23db 

where ( I34bb is derived by using the Schur complement through 
the auxiliary variable Y. 

The role of ^rHQfe — QaT Hi? i s to render the cost 
in ( I34ab strictly convex and coercive. Moreover, for small 
values of T k , the optimization variable Q k is forced to stay 
"close" to Qi" 13 obtained at the previous iteration, thereby 
improving the stability of the iterates [34 Ch. 6]. Centralized 
and distributed schemes with the proximal point regularization 
are given by Algorithms Q] and |2j respectively, with problem 
(P6) replacing (P5). Convergence of the resulting schemes is 
established in the following theorem. To avoid ambiguity, these 
proximal point-based algorithms will be hereafter referred as 
Algorithms Q~fP) and^P), respectively. 

Theorem 2. Suppose that the sequence {Q*-™- 1 } generated by 
Algorithm 1(P) (Algorithm 2(P)) has a limit point. Then, every 
limit point is a stationary point of (P3). 

Proof: The Gauss-Seidel method with a proximal point 
regularization converges without any underlying convexity 
assumptions 1351 . A modified version of the proof is reported 
here, where the local convex approximation ( TT3l and the 
peculiarities of the problem at hand are leveraged to establish 
not only convergence of the algorithm, but also optimality of 
the obtained solution. 
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Assume there exists a subsequence {Q^ n ^\j — 1,2,...} 
converging to a limit point Q := (Qi, . . . , Qk)- Let Qfc +1 ' > 
be obtained as 

n («+i) — 

*«4fc ■— 

argmax U k ( Q fc ; Q(" +1) , . . . , Q<£?\ qft, . . . , 



QkeQ k 



-^IIQ fe -Qi n) lll- (35) 

Thus, it follows that [cf. ([26t l 

UiiQ^-A^) > ^(Q^; QL" j) ) 

+ J_||Q^+i)_Q^)|||. (36) 

Going along the lines of the proof of Theorem Q] it holds that 



lim ^(Q^+'^Q^ ) = lim Z4(Qp' J ! Q-T ) = ^i(Q). 

(37) 



Therefore, taking the limit as j — > oo in (136t . one arrives at 



lim IIQ^-Q^^lll 

j-s-oo 







(38) 



which implies that Q^™ 3+1 ^ also converges to Qi. 

Since Q^™ 3+1 ^ is generated as in d35l >, it satisfies the 
optimality condition 

m Tr 



V!^(Q^ +1) ; Q^) - -(Qi l3+1) ~ Qi^) 
n 



n 



(Qi - Qi 



< 0, V Qi g Ql (39) 



Taking the limit as j — !> oo in d39l , and using again the fact 
that ViWi(Q) = ViWi(Q), we obtain 

5R{Tr{Vi^(Q) w (Qi -Qi)}} < 0, V Qi 6 Qi . (40) 

Then, repeating the same argument for all k S K, leads to 

5i{Tr{V fc W(Q) w (Q fe -Q fc )}} < 0, V Q fc e Q fc (41) 

which shows that the limit point Q is also a stationary point. 

■ 

As asserted in Theorem Algorithms 1(P) and 2(P) con- 
verge to a stationary point of (P3) for any possible antenna 
configuration. The price to pay however, is a possibly slower 
convergence rate that is common to proximal point-based 
methods [34, Ch. 6] (see also the numerical tests in SectionlVT). 
For this reason, the proximal point-based method should be 
used in either a centralized or a distributed setup whenever 
the number of transmit-antennas exceeds that of receive- 
antennas in at least one transmitter-receiver pair. In this case, 
Algorithms 1(P) and 2(P) ensure first-order optimality of the 
solution obtained. When < Nk, for all k £ K, the two 
solvers have complementary strengths in convergence rate and 
computational complexity. Specifically, Algorithms 1 and 2 
require the rank of all CR direct channel matrices {Hfc^} 
beforehand, which can be computationally burdensome, es- 
pecially for a high number of antenna elements. If the rank 
determination can be afforded, and the convergence rate is at 
a premium, then Algorithms 1 and 2 should be utilized. 



IV. Aggregate Interference Constraints 



Suppose now that the individual interference budgets {i™ x } 
are not available a priori. Then, the aggregate interference 
power {i max } has to be divided among transmit-CRs by the 
resource allocation scheme in order for the overall system 
performance to be optimized. Accordingly, (P3) is modified 
as follows to incorporate a robust constraint on the total 
interference power inflicted to the PU node: 



A' 



(P7) max y>fe({Q fe }) 



s. t. Tr{Q fe } < pT\ k G K 



(42a) 
(42b) 



K 



]TTr{G fc Q fc G«} < i max , V AG, e C7 fc , k e AC. (42c) 



fe=i 



The new interference constraint ( I42ct couples the CR nodes 
(or, more precisely, the subset of transmit-CR nodes in the 
proximity of the PU receiver). Thus, the overhead of message 
passing increases since cooperation among coupled CR nodes 
is needed. 

A common technique for dealing with coupled constraints is 
the dual decomposition method l29l . which facilitates evalua- 
tion of the dual function by dualizing the coupled constraints. 
However, since (P7) is non-convex and non-separable, the 
duality gap is generally non-zero. Thus, the primal variables 
obtained during the intermediate iterates may not be feasible, 
i.e., transmit-covariances can possibly lead to violation of the 
interference constraint. Since the ultimate goal is to design an 
on-line algorithm where ( I42cb must be satisfied during network 
operation, the primal decomposition technique is well moti- 
vated to cope with the coupled interference constraints iflOll . To 
this end, consider introducing two sets of auxiliary variables 
{tfc} and {tk} in problem (P7), which is equivalently re- 
formulated as 



(P8) 



A 

max VTr{H fc , fc Q fc H« fc (H fe , fe Q fe H« fc + R fc , fc ) _1 } 

(43a) 

s. t. Tr{Q fc } < p^ x , keK, (43b) 
Tr{G fe Q fe G«} < t k , V AG fc e Q k , k e K (43c) 
< t k < L k , keJC (43d) 

A 

J2 * ^ tmax • ( 43e > 
/c=i 

For fixed {i k }, the inner maximization subproblem turns out 
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to be 
(P9) p({t fc }):= 

K 

max VTr{H M Q t H« fc (H fc , fe Q fc H* fc + R fe , fe ) _1 } 

(44a) 

s. t. Tr{Q fc } < * G £ (44b) 

Tr{G fe Q fe G«} < t fc , V AG fc G fc G IC (44c) 
< t fc < i k , keJC (44d) 

which, as discussed in preceding sections, can be solved using 
the block coordinate ascent algorithm (or its proximal point 
version) in either a centralized or a distributed fashion. After 
solving (P9) for a given set {ik}, the per-CR interference 
budgets {tfc} are updated by the following master problem: 

(P10) max p({t k }) (45a) 
s. t. {i k } G 1 (45b) 
with the simplex set I given by 

X := jujkfe > 0,J2^k < i max | • (46) 

Overall, the primal decomposition method solves (P8) by 
iteratively solving (P9) and (P10). Notice that the master prob- 
lem (P10) dynamically divides the total interference budget 
t mal among CR transmitters, so as to find the best allocation 
of resources that maximizes the overall system performance. 
Using the block coordinate ascent algorithm, the fc-th transmit- 
covariance matrix is obtained by solving the following 
problem [cf. Algorithm 2] 

CPU) Pk(tk):= 

max u fc (Q fc ,Q_ fe )+Tr{D«Q fc } (47a) 

s. t. Tr{Q fe }< P r x (47b) 
Tr{G fc Q fe G^}<i fc ,VAG fc eg fc (47c) 
t k < t* (47d) 

where the proximal point-based regularization term is added if 
Algorithm 2(P) is implemented. Since (PI 1) is a convex prob- 
lem, it can be seen that the subgradient of pu^k) with respect 
to ik is the optimal Lagrange multiplier A& corresponding to 
the constraint ( I47db [29, Chap. 5]. Thus, it becomes possible 
to utilize the subgradient projection method to solve the 
master problem. Strictly speaking, due to the non-convexity 
of the original objective ( |43at , primal decomposition method 
leveraging the subgradient algorithm is not an exact, but rather 
an approximate (and simple) approach to solve (P8). However, 
because (147 at is a tight concave lower bound of (143 al l around 
the approximating feasible point, p({i k }) is well-approximated 
by Pfc(tfc) as {Qt } approaches the optimal value {Q^, pt }. 
Hence, A^ also comes "very close" to the true subgradient of 
Pfe(j>fe}) with respect to ik- Therefore, at iteration £ of the 



Algorithm 3 Distributed on-line robust sum-MSE minimiza- 
tion with aggregate interference constraint 

l: Initialize Q[ 0) (0) = 0, and i k (0) = L mm /K, V k G JC 

2: repeat (£ = 1,2,...) 

3: [CRs]: Solve (P9) via Algorithm 2 [Algorithm 2(P)] 
4: [CRs]: Transmit {Xk{£)} to the cluster-head node 
5: [Cluster-head node]: Update {tk{£ + 1)} via (l48l 
6: [CRs]: Receive 1)} from the cluster-head node 

7: until u k (Q (n) {£)) - u k (q^(1 - 1)) < v, V k G K 



primal decomposition method, the subgradient projection up- 
dating the interference budgets t := [ti, ta> • • • > lk] T becomes 

l(£ + 1) = Proj I [t(£) + s(£)\(£)} (48) 

where A := [Ai, A2, . . . , \k] T \ s(£) is a positive step size; 
Projj[ ] denotes projection onto the convex feasible set I. 
Projection onto the simplex set in d46b is a computationally- 
affordable operation that can be efficiently implemented as in 
e.g., EH. 

Once (P9) is solved distributedly, each CR that is coupled 
by the interference constraint has to transmit the local scalar 
Lagrange multiplier \k{£) to a cluster-head CR node. This 
node, in turn, will update {tk{£ + 1)} and will feed these 
quantities back to the CRs. The resulting on-line distributed 
scheme is tabulated as Algorithm 3. Notice however that in 
order for the overall algorithm to adapt to possibly slowly 
varying channels, operation (l48l i can be computed at the end 
of each cycle of the block coordinate ascent algorithm, rather 
than wait for its convergence. 

V. Simulations 

In this section, numerical tests are performed to verify the 
performance merits of the novel design. The path loss obeys 
the model d~ n , with d the distance between nodes, and rj = 
3.5. A flat Rayleigh fading model is employed. For simplicity, 
the distances of links XJ\ — > UJ1 are all set to d k ,k = 30 
m; for the interfering links {Ul — > UJ,j ^ k) distances 
are uniformly distributed over the interval 30 — 100 m. As 
for the distances between CR transmitters and PU receivers, 
two different cases are considered: (cl) the PU receivers are 
located at a distance from the CRs that is uniformly distributed 
over 70 — 100 m; and, (c2) the CR-to-PU distances are 
uniformly distributed over 30 — 100 m. Finally, the maximum 
transmit-power and the noise power are identical for all CRs. 
For the proximal point-based algorithm, the penalty factors 
{rk} are selected equal to 0.1. 

To validate the effect of the robust interference constraint, 
the cumulative distribution functions (CDF) of the interference 
power at the PU are depicted in Fig. [2] Four CR pairs and 
one PU receiver are considered, all equipped with 2 antennas. 
The maximum transmit-powers and noise powers are set so 
that the (maximum) signal-to-noise ratio (SNR) defined as 
SNR := pr x (dkl)/^l e q uals 15 dB - The total interference 
threshold is set to t max = 4 • 10~ 7 W and, for simplicity, 
it is equally split among the CR transmitters. The channel 
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uncertainty is set to e| = p ■ \\G k \\ 2 F 02), with p = 0.05. 
CDF curves are obtained using 2,000 Monte Carlo runs. In 
each run, independent channel realizations are generated. The 
Matlab-based package CVX [37 1 along with SeDuMi J38) are 
used to solve the proposed robust beamforming problems. 

The trajectories provided in Fig. [2] refer to the block co- 
ordinate ascent (BCA) algorithm described in Section [Till the 
one with the proximal point-based regularization (proximal- 
BCA) explained in Section IIII-CI and the non-robust solver 
of (P2), where the estimates {Gk} are used in place of 
the true channels {G&}. Furthermore, the green trajectory 
corresponds to (P8), where the subgradient projection (|48T > is 
implemented at the end of each BCA cycle, which includes K 
updates of Qfc for k = 1, . . . ,K. As expected, the proposed 
robust schemes enforce the interference constraint strictly in 
both scenarios (cl) and (c2). In fact, the interference never 
exceeds the tolerable limit shown as the vertical red solid 
line in Fig. [2] The CDFs corresponding to the proposed 
BCA and its proximal counterpart nearly coincide. In fact, 
the two algorithms frequently converge to identical stationary 
points in this particular simulation setup. Notice that with the 
primal decomposition approach the beamforming strategy is 
less conservative. On the contrary, the non-robust approach 
frequently violates the interference limit (more than 30% of 
the time). Finally, comparing Fig. |2(a)| with Fig. |2(b)| one 
notices that the interference inflicted to the PU under (cl) and 
the one under (c2) are approximately of the same order. Since 
in the second case the CR-to-PU distances are smaller, the 



CR transmitters lower their transmit-powers to protect the PU 
robustly. 

Convergence of the proposed algorithms with given channel 
realizations and over variable SNRs is illustrated in Fig. [3] It 
is clearly seen that the total MSEs decrease monotonically 
across fast-converging iterations, and speed is roughly iden- 
tical in (cl) and (c2). As expected, the proximal point-based 
algorithm exhibits a slightly slower convergence rate. Notice 
also that the primal decomposition method returns improved 
operational points, especially for medium and low SNR values. 
Furthermore, the gap between the sum-MSEs obtained with 
and without the primal decomposition scheme is more evident 
under (c2). Clearly, the sum-MSEs at convergence in (c2) are 
higher than the counterparts of (cl). This is because CRs are 
constrained to use a relatively lower transmit-power in order to 
enforce the robust interference constraints; this, in turn, leads 
to higher sum-MSEs and may reduce the quality of the CR- 
to-CR communications. 

In Fig.|4j the achieved sum-MSE at convergence is reported 
as a function of the total interference threshold. Two sizes 
of the uncertainty region are considered with p = 0.05 and 
p = 0.1. Focusing on the first case, it can be seen that the 
two achieved sum-MSEs first monotonically decrease as the 
interference threshold increases, and subsequently they remain 
approximately constant. Specifically, for smaller t max , the 
transmit-CRs are confined to relatively low transmit-powers 
in order to satisfy the interference constraint. On the other 
hand, for high values of i max , the interference constraint is no 
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«— BCA, p = 0.1 

BCA w/ primal decomp., p = 0.1 
a— BCA, p = 0.05 

BCA w/ primal decomp., p = 0.05 
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Fig. 4. Achieved sum-MSE as a function of i ma \ for SNR = 10 dB. 
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Fig. 6. Achieved sum-MSE for SNR = 15 dB. 
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Fig. 5. Achieved sum-MSE for SNR = 15 dB. 



longer a concern, and the attainable sum-MSEs are mainly due 
to CR self-interference. Notice also that for p = 0.1 the sum- 
MSEs are clearly higher, although they present a trend similar 
to the previous case. This is because the uncertainty region 
in ( 1 1 2cb becomes larger, which results in a higher sum-MSE. 

In order to compare performance of the proposed algo- 
rithms, the total MSE obtained at convergence is depicted 
in Fig. [5] for 50 different experiments. In each experiment, 
independent channel realizations are generated. The SNR is set 
to 15 dB. It is clearly seen that the objectives values of the two 
proposed methods often coincide. The differences presented in 
a few experiments are caused by convergence to two different 
stationary points. In this case, it is certainly convenient to 
employ the first algorithm, as it ensures faster convergence 
(see Fig. |3(a)[ ) without appreciable variations in the overall 
MSE. Notice that a smaller mean-square error can be obtained 
by resorting to the primal decomposition technique. 

In Fig. |6l the simulation setup involves 8 CR pairs and one 
PU receiver. The CR transmitters have 4 antennas, while the 
receiving CRs and the PU are equipped with 2 antennas. The 
distances dk,k are set to 50 m, while {dkj}k^j distances are 
uniformly distributed in the interval between 30 and 250 m. 
Finally, CR-to-PU distances are uniformly distributed between 
100 and 200 m. Clearly, matrices {H/^,} here do not have 
full column rank. It is observed that about 10% of the times 
the proximal point based algorithm yields smaller values 
of the sum-MSE than Algorithm 1. This demonstrates that 
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Proximal-BCA, setup of Fig. 5 

BCA w/ primal dec, setup of Fig. 5 
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Proximal-BCA, setup of Fig. 6 

BCA w/ primal dec, setup of Fig. 6 



Fig. 7. CDF of sum-MSE gaps (relative to the BCA) using proximal-BCA 
(blue) and primal decomposition (red): MSE(BCA) — MSE(proximal-BCA) 
and MSE(BCA) - MSE(BCA-primal decomp.). 



Algorithm 1 may not converge to a stationary point, or, it 
returns an MSE that is likely to be worse than that of the 
proximal point-based scheme. 

Fig- depicts the CDFs of the difference between the sum- 
MSE obtained with BCA, along with the ones obtained with 
proximal-BCA and with the primal decomposition method. 
The simulation setups of Figs. [5] and [6] are considered. In the 
first case, it can be seen that for over 80% of the trials the BCA 
and proximal-BCA methods yield exactly the same solution. 
Moreover, BCA with primal decomposition performs better 
than the BCA method about 90% of the time. Specifically, the 
gain can be up to 0.765, which corresponds to approximately 
14% of the average sum MSE of the BCA. In the second case, 
the proximal-BCA returns a smaller sum-MSE with higher 
frequency. 

VI. Concluding summary 

Two beamforming schemes were introduced for underlay 
MIMO CR systems in the presence of uncertain CR-to-PU 
propagation channels. Robust interference constraints were 
derived by employing a norm-bounded channel uncertainty 
model, which captures errors in the channel estimation phase, 
or, random fading effects around the deterministic path loss. 
Accordingly, a robust beamforming design approach was for- 
mulated to minimize the total MSE in the information symbol 
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reconstruction, while ensuring protection of the primary sys- 
tem. In order to solve the formulated non-convex optimization 
problem, a cyclic block coordinate ascent algorithm was devel- 
oped, and its convergence to a stationary point was established 
when all CR-to-CR direct channel matrices have full column 
rank. A second algorithm based on a proximal point regulariza- 
tion technique was also developed. Although slower than the 
first, the proximal point-based scheme was shown capable of 
converging to a stationary point even for rank-deficient channel 
matrices. The two solutions offer complementary strengths 
as far as convergence rate, computational complexity, and 
MSE optimality are concerned. They can both afford on-line 
distributed implementations. Finally, a primal decomposition 
technique was employed to approximately solve the robust 
beamforming problem with coupled interference constraints. 
The developed centralized and distributed algorithms are also 
suitable for non-CR MIMO ad-hoc networks as well as for 
conventional downlink or uplink multi-antenna cellular sys- 
tems. 



Appendix 

A. Derivation of the complex gradient matrix d!61 > 

Using that — — - = A 113911 , and letting [A] mn denote 
the (m, n)-th entry of matrix A, it follows that 

9[B«] st _aTr{e«H J , fc Q«H«e t } 



dQl 



d[Bf] 
d[Q 



- H ,-.fe e t e s H i,fc 



(49) 



e m H j,fc e * e s tlj^en = e s Hj, k e n e m H j k et 



(50) 



which can be written in a compact form as 

3[Q*f m „ = H j,fc e " e m H jV Then > the identity 

J£r = -X- 1 ( a( x-i). ) x_1 G3> which holds for 
any Hermitian positive definite matrix X, is used to obtain 



dB* 



(51) 



Using now the chain rule, one arrives at 

T 



Tr 



OB? 



d[Q* k } m n " \ \dBf J d[Ql iW! , 

= Tr{ e*H&B: V ,B (52) 

which readily leads to the desired result 

du 3 _ ttW R-i^.n-ii 



dQt 



(53) 



B. Proof of Lemma \2\ 

First, convexity of Q k can be readily proved by the def- 
inition of a convex set ifTSl Ch. 2]. Re- write the function 
«i(Qfc,Q-fc) as [cf. ©, d] 

Uj(Q k , Q_ fe ) = Tr{vy 2 PTi(Q fe )V 3 1/2 } (54) 

where 

Pj(Qfc) = TT, ;.Q/,1I J', + If , ,Q,nj', + rrjT.y, 

i^k 

is an affine map with respect to Q^. Since Uj is convex 
in Pj 11401 Theorem 2], and convexity is preserved under 
affine mappings and nonnegative weighted-sums lfl8l Ch. 3], 
it follows that /fc(Qfe,Q-fc) is convex in Qfe. 

C. Proof of Lemma \3\ 

First, notice that the objective function d!5ab can be re- 
written as 

Uk(Qk) =A^+Tr{D«Q fc } 

-Tr{R fc , fe (H fe , fe Q fc H« fc +R M ) _1 } . (55) 

Then, it suffices to prove strict convexity in Qk of the third 
term on the right hand side of d55l l. This is equivalent to 
showing that (subscripts are dropped for brevity) 



J(t) := Tr |r (HQH w + R) H 



(56) 



is strictly convex in t G {f |Q := X + iY g Q} for any given 
X e H" x " and nonzero Y e H nx ™. 

To this end, consider the second-order derivative of J(i), 
which is given by 

J(t) = 2Tr{CRCLCL} (57) 

where C := (R + H(X + tY)H' H ) 1 and L := HYH M . 
Note that matrix CRC is Hermitian positive definite, since 
C and R are Hermitian positive definite too. With H full 
column rank, it readily follows that L ^ for any Y / 0. 
This ensures that the Hermitian positive semi-definite matrix 
LCL is not an all-zero matrix, i.e., LCL ^ 0. 

Let v\ > V2 > ■ ■ ■ > vn > and /xj > /Z2 > ■ ■ ■ > 
Hn > denote the eigenvalues of matrices CRC and LCL, 
respectively. Since matrix LCL ^ 0, /Ltj is strictly positive, 
and thus 



N 
i=l 

> 2vn[1i > 



(58a) 
(58b) 



where (15 8 at follows from von Neumann's trace inequal- 
ity BP . Finally, (I58bb shows the strong convexity (and hence 
strict convexity) of J(t). 

For completeness, we provide an alternative proof of 
the lemma. With some manipulations, function h(Q) := 
Tr |r (HQH w + R) 1 j can be re-expressed as 

h(Q) = ^R-^HQH^R- 1 / 2 ) 



Tr {(i + R-^HQH^R- 1 / 2 ) 



(59) 
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where g(X) := Tr {(I + X) _1 |. Let Ai(X), . . . , A„(X) de- 
note again the eigenvalues of a matrix X. Note that the spectral 
function g(X) = s(A(X)) := J2i (l+X~(X)) * s strictly convex 
if and only if the corresponding symmetric function s(-) is 
strictly convex ll42ll . To this end, the strict convexity of 
for x > implies the strict convexity of s(-), and thus of <?(X). 
Under the condition of full column rank of H, we will show 
that strict convexity is preserved under the linear mapping 
in ( |59l l. Specifically, define 

Q 4 := R - 1/2 HQiH w R _1/2 , i = 1,2. 
Then, for any Qi ^ Q 2 G <2 and < A < 1, we have that 
/i(AQi + (1 - A)Q 2 ) = ,g(AQi + (1 - A)Q 2 ) (60a) 
< A. 9 (Qi) + (1 - X)g(Q 2 ) (60b) 
= Xh(Qx) + (1 - A)/i(Q 2 ) (60c) 

where (I60bb follows from the strict convexity of g(-), and the 
fact that Qi 7^ Q 2 holds for any Qi 7^ Q 2 , since H is full 
column rank. 
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